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We study the Stell-Hemmer potential using both analytic (exact Id and approximate 2d) solutions 
and numerical 2d simulations. We observe in the liquid phase an anomalous decrease in specific 
volume and isothermal compressibility upon heating, and an anomalous increase in the diffusion 
coefficient with pressure. We relate the anomalies to the existence of two different local structures 
in the liquid phase. Our results are consistent with the possibility of a low temperature/high pressure 
liquid-liquid phase transition. 
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In their pioneering work, Stell and Hemmer proposed 
the possibility of a new critical point in addition to the 
normal liquid-gas critical point for potentials that have 
a region of negative curvature in their repulsive core 
(henceforth referred to as core-softened potentials) [Q. 
They also pointed out that for the \d model with a long 
long range attractive tail, the isobaric thermal expan- 
sion coefficient, ap = {dV/dT)p (where V,T and P 
are the volume, temperature and pressure) can take an 
anomalous negative value. Debenedetti et al., using ther- 
modynamic arguments, pointed out that the existence of 
a "softened core" can lead to ap < ||^. 

Here we further investigate properties of core-softened 
potential fluids. We first study the properties of the ID 



fluid using an exact solution. We then investigate the 
behavior of the 2D fluid, initially by an approximate 
solution provided by cell theory method and finally by 
performing molecular dynamics simulation of the fluid. 
The discrete form of the potential that we study is 



u{r) — 



oo < r < a 

— Ae a < r < b 

— e b < r < c 

r > c 



(1) 



with r being the inter-particle distance and A < 1 
(Fig. 0(a)) [||. The model is exactly solvable in Id, fol- 
lowing the methods of [^0, and the equation of state 
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Here V/N = ^ is the average distance between near- 
est neighbors, H^ = e^^^^ {x — a,b,c), W = e^'^ and 
/3 = ksT. The isobars (Fig. |l|(b)) exhibit two differ- 
ent types of behavior. For all P larger or equal to an 
upper boundary pressure Pup, £ = a aX T = 0, and £ 
increases monotonically with T. For P < P^p, £ = b a.t 
T = 0. The isobars show a maximum and a minimum in 
£, which correspond respectively to points of minimum 
and maximum density]^, bounding a density anomaly 
{ap < 0) region There is a discontinuity in £ at 

P = Pup along the T = isotherm. 

Next we study the isothermal compressibility Kt = 
-y-i {dV/dP)j,j^. We use Eq.(|) to calculate Kt along 
isobars (Fig. |l|(c)). The graphs show an anomalous re- 
gion in which Kt decreases upon heating (for simple liq- 
uids Kt increases with T). We find the maximum value 



of Kt grows as P increases towards Pup, and Kt di- 
verges as 1 /T when we approach the point C" with coor- 
dinates (T — 0,P = Pup) which we interpret as a critical 
point pil. Further, the locus of Kt extrema joins the 
point C'(Fig. 0(d)). 

We also study the Tmd locus (Fig. 0(d)) and note that 
the locus of Kt extrema intersects the Tf^ip) locus at its 
infinite slope point, a result that is thermodynamically 
required JT2|. Such a point on the Tmd has been ob- 
served in simulations which support the existence of a 
liquid- liquid phase transition in supercooled water . 

We next consider the 2d case, for which an exact solu- 
tion does not exist. We use the spherical Lennard- Jones 
and Devonshire cell theory method which assumes 
that each particle is confined to a circular cell, whose 
radius is determined by the average area per particle, v. 
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This method neglects the correlation between the posi- 
tions of different particles and assumes that the poten- 
tial acting on each particle is a result of interacting with 
all its nearest neighbors smeared around its cell. The 
Helmholtz free energy per cell is 

h{v,T) - h,iv,T) - kBT\n{vf{v,T)/v), (3) 

where hi{v, T) is the ideal gas free energy and Vf{v,T) is 
the free volume defined as 

Vf{v,T)^ [ e-^"(-)dx (4) 
J cell 

with the core-softened potential used for u(x) . For each 
value (P, T), we find the value of v{P, T) by minimizing 
h{v,T) — Pv. The resulting phase diagram (Fig. |^) has 
two lines of first order phase transition, a low pressure 
line which is the liquid-gas phase transition line termi- 
nating at a critical point C, and a high pressure line that 
separates a low-density liquid (LDL) and a high-density 
liquid (HDL) and terminates in a critical point C". This 
picture holds both for the discrete and smooth versions 
of the potential. We note that the presence of C" would 
imply an anomalous increase in Kt upon cooling when 
C is approached from higher temperatures. 

In order to further investigate the system in the liquid 
phase we rely on numerical molecular dynamics (MD) 
study of the discrete and smooth versions of the poten- 
tial (Fig. [^(a)). We perform 2d MD simulations for a 
system composed of N circular disks, in a rectangular 
box of area A. To each disk we assign a radius equal to 
half the hard core diameter a, and define the density p as 
the ratio of the total area of all the disks to the area of 
the box. For the discrete version of the potential we use 
the collision table technique jl6j for N — 896 disks and 
for the smooth version of the potential, we use the ve- 
locity Verlet integrator method jl^] with N — 2500. We 
set the mass of the particles to be unity, while the units 
of length and energy are scaled by e and a respectively. 
We choose the time step 6t = 0.01 ||r^. For each T, we 
first slowly thermalize the system, using the Berendsen 
method of rescaling the kinetic energy jl^, after which 
we perform the simulation at constant N, A and energy 
E. We fix p by fixing A and we start from T = 1, low- 
ering T down to 0.4 (in steps of AT = 0.1 for larger T 
and AT = 0.05 and 0.025 in the vicinity of freezing). As 
the initial configuration for each (p, T) , we choose the 
equilibrated configuration of (p, T -1- AT). We simulate 
state points along constant p paths (isochores) (Fig. ||(a)) 
and also along constant P paths (isobars) We use 

the isobar results to check the values of p{T, P) calcu- 
lated from isochores, as well as to find Kt along isobars 
(Fig. 1(b)) (20|. 

The MD results are qualitatively equivalent for the dis- 
crete and smooth versions of the potential. For the re- 
sults of Fig. I we have used the smooth version. The 
minima in the P versus T isochores (Fig. |(a)) corre- 
spond to density maximum points We also find 



that along some of the isobars, Kt increases upon cool- 
ing (Fig. g(b)) H^. The Tmd locus possesses a point 
with infinite slope (as in Id), which we verify by find- 
ing the intersection of the locus of Kt minima with the 
Tmd line in Fig. ||(a)||l^. If we assume that a metastable 
liquid critical point exists below freezing, then by fitting 
the Kt graph to a power law divergence, we can estimate 
the critical point to be in the region 0.3 < T < 0.5 and 
1.0 < P < 1.5 which is in agreement with the cell-theory 
approximation. 

We also study the effect of pressure on diffusion, and 
find that along some isotherms, increasing pressure in- 
creases D, while for simple liquids increasing pressure 
decreases D (Fig. ^(c)). This anomaly occurs in the same 
region of phase diagram where the density and isothermal 
compressibility anomaly is observed. 

The anomalies can be related to the interplay between 
two local structures, an open structure in which the near- 
est neighbor particles are typically at a distance b and a 
denser structure in which the nearest neighbors penetrate 
into the softened core and are typically at a smaller dis- 
tance a. The configurations are determined by the min- 
ima of the Gibbs free energy, G{T, P) = U + PV -TS 
(where G, U and S are the Gibbs free energy, internal 
energy, and entropy). Fig. |^(a) shows the Id free energy 
at T = for two different values of P. The qualita- 
tive shape should not change for higher dimensions and 
small T. For low pressures at small T, the open struc- 
ture is favored by the free energy. Increasing T for these 
pressures will increase local fluctuations in the form of 
dense structures which can lead to an overall contraction 
of the system upon heating, causing a density anomaly. 
Increasing P raises the relative free energy of the open 
structure, until the dense structure will be the favored lo- 
cal structure, as seen for P > P„p in Fig. ^(a). At small 
T this can lead to a first order pressure driven transi- 
tion ("core-collapse" Q]), while for large T the transition 
is continuous. At T = 0, the value of Pup where the 
transition occurs is where G(T = 0, P) is the same for 
the open and dense structures, so 

T-j Uopen Udense /^n 

open ^ dense 

From Eq.(0) and Eq.(|) we find P„p = (1 - A)e/(6 - a) 
in Id, which can also be derived from Eq.(|[)||lJ]. 

To examine the transition from the open structure to 
dense structure in 2d, we study the pair distribution func- 
tion g{r) for the MD configurations (Fig. §(b)). The first 
peak in g{r) splits into two subpeaks, which correspond 
to the locations of the nearest neighbors in the dense 
and open structures. As T increases, the open structure 
subpeak decreases while the dense structure subpeak in- 
creases. We observe the same change with P along the 
liquid isotherms for small T. The uniform value of g{r) 
for large r confirms that all the state points of Fig. ||(a) 
are in the liquid state. 
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FIG. 1. (a) General form for the core-softened potential 
studied here. The length parameters a, b, c and energy param- 
eters e, A are shown. The dashed curve is the smooth version 
[^. (b) Isobars (the average distance per particle, £ versus T 
at constant P) for the discrete Id core-softened potential|^] 
with Pup = 2.5 in agreement with Eq.(|^. The Tmd point 
is marked by an open circle, (c) Isothermal compressibility 
for the discrete potential along different isobars, with their 
maxima marked by filled circles. The Kt for Pup diverges as 
1/T. (d) The loci of Tmd and Kt extrema for the discrete 
potential. 
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FIG. 2. (a) Solid lines correspond to the phase diagram 
from the cell theory method for the 2d discrete potential with 
the same parameters as in Fig. |l| (solid curve inset) . The lower 
solid line corresponds to the liquid-gas (L-G) phase transition, 
with a critical point outside the range of the graph. The upper 
solid line corresponds to the HDL to LDL phase transition, 
the slope of which changes from negative to positive. The 
Clausius-Clapeyron equation relates the slope of a first or- 
der transition line to the difference between the entropy and 
volume of the phases by dP/dT = AS/AV. In the case of wa- 
ter, the open LDL structure produced by highly directional 
hydrogen bonding has a lower entropy. We can mimic this sit- 
uation for our radially symmetric potential by narrowing the 
well which greatly reduces the positively sloped part of the 
transition line; the dashed line corresponds to the HDL-LDL 
transition for a discrete potential with c = 1.40001. (b) Simi- 
lar phase diagrams for the smooth version of the potential in 
Fig. |l|(a) (solid lines) . The dashed line correspond to a poten- 
tial with a much narrower well (w — 10, 000). 
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FIG. 3. MD results for the smooth potential with the 
same parameters as in Fig. |l||P. (a) Constant density curves 
with, from bottom to top, densities between 0.46 to 0.56 in 
steps of 0.1. The open circles mark Tmd and the dashed line 
is the locus of Kt minima from part (b). The thick gray line 
is the approximate locus of the freezing points, (b) Isothermal 
compressibility along isobars. Except for the P = 0.25 isobar, 
the graphs show anomalous decrease upon heating, (c) Diffu- 
sion coefficient D (slope of the mean square displacement as a 
function of time) for different pressures at T = 0.65, showing 
an anomalous increase in the P < 1.5 range. For comparison, 
we show D/4 for a Lennard- Jones liquid at T = 0.7, from 
simulation of 2304 disks. The high pressure zero values in 
both graphs correspond to a solid phase. 
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FIG. 4. (a) The Id Gibbs free energy at T = as a func- 
tion of the extra "degree of freedom" £, for the discrete form 
of potential in Fig. [lj(a). The equilibrium value of 1{P) is 
determined as the absolute minimum of this function, which 
is located at ^ = fe below Pup and a.t £ = a above Pup- (b) 
Pair distribution function for T = 0.6 and T — 1.0 on the 
p = 0.48 isochore of Fig. ^(a). The uniform g{r) for large r 
is a sign of the liquid phase. The inset is a blow up of the 
first (split) peak. The thick solid curve is for T = 1.0 and 
the dashed curve is for T — 0.6, indicating that increasing 
T lowers the total number of particles in the open structure 
(r « 1.5) and increases the number of particles in the dense 
structure (r ~ 1.1). 



